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АВ5ТКАСТ 


A variety of simple analytical models for increasing, 
decreasing and "bath tub"-type failure rates are discussed. 
The purpose of this thesis is to develop analytical hazard 
representations for use in reliability and maintainability 
studies, and to evaluate them in use for data analysis. 
Verification of the model was accomplished by computer simu- 
lation. They were applied to human mortality and other 


failure time data. 
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I. INTRODUCTION AND SUMMARY 


A. INTRODUCTION 

The failure rate function, or hazard function (hazard 
for short) may be described as the conditional probability 
of an equipment's failing at operating age t; having sur- 
vived to that age. The reliabilities of a variety of elec- 
tronic and mechanical items are conveniently and naturally 
described in terms of the appropriate hazard function, and 
—7 7:7. longevity of human beings. The term force of 
Mee e replaces hazard in the latter context. 

This paper is devoted to a study of several simple 
analytical representations for hazard functions. These 
representations are in turn based upon representations of 
random variables having certain required properties, in 
terms of others having familiar distributions--in particular 
the exponential. Similar ideas are due to Tukey [1], and 
recently have been examined by Parzen [2]. The hazard 
representations proposed are quite expeditiously used in 
simulation studies, e.g. of system reliability or avail- 
ability in terms of component lifetimes. They may also be 
used in data analysis studies, in order to parsimoniousiy 
describe data sets in terms of perturbations of convenient 
and familiar standard distributions. Their use in data 
analysis and simulation is also described in Gaver, Laven- 


berg, and Price [3], and in Gaver and Chu [4]. 








B. SYSTEM FAILURE PATTERNS 
It is plausible to think that the time series of fail- 


ures in a system may involve these stages. 


Early failures. There may be a relatively large number 
of failures soon after a system is introduced because of 
design defects, production errors, or errors stemming from 
maintenance personnel inexperience. This situation is 
characterized by a hazard function that is initially large, 
but that decreases with time. "Infant mortality" is in 


evidence. 


Random Failures. Following the early failure period 
there may be a period during which failures occur at an 
essentially constant rate for a rather prolonged time. 
During this period the hazard function is nearly constant, 
So the times between failures are close to being exponen- 
tially distributed. The effect of age or wearout is not 


Mer apparent. 


Wearout Failures. Eventually following the period 
Suringswhich a Constant hazard is evident there is likely 
to be a period of ever-increasing failure rate caused by 


wearout of system components. 


A graphical representation of a hazard function that 
exhibits the behavior described is given below. Note that 


it has the legendary "bath tub" shape. 
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Some comments on the above follow: 

The term "failure" may refer to an event that is 
analogous to human death, after which the entire system 
is replaced. On the other hand repair or component replace- 
ment may occur after failure: the system is only repaired, 
not entirely EECH In the former case, a hazard func- 
tion of the kind depicted in Figure l applies to each system 
event ("death"); when the system is installed (or is born), 
that hazard operates starting from scratch at t = 0 until 
system failure (human death, for instance), after which a 
similar hazard goes into effect, starting once again from 
zero. In the latter case, in which repair of a component 
occurs, a hazard function like that of Figure l applies at 


mb 0 Dub after the first event ("failure") at t a re- 


1 
pair action is accomplished. The same hazard operates for 
E» t. until the next event at 6. ә Еј, and so on.  Inter- 
mediate situations may be envisioned, in which after event 

n at t, the hazard governing system failure п+1 starts 


meme to -— t, 0 < t. «X t. 
n n n n 
Although there is reason to assume that hazards some- 
what like that of Figure l occur in general for systems, 
the possibility exists that the system hazard is "bumpy" 
because wearout failures of components or subsystems may 
well occur at intermediate times. 
If the theory is applied to systems with little or no 


wearout propensity, as should be the case when dealing with 


computer software modules, then the hazard function may 
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well exhibit the initial falloff of Figure 1 but not the 


rise at later times. In fact, a constant decline as bugs 


are found and removed could be (optimistically) anticipated 


for software. The right-hand side of the bath tub vanishes, 


and the picture is that of a ski slope. 


El 





II. ANALYTICAL HAZARD REPRESENTATIONS 


A. MODELS FOR THE HAZARD FUNCTION 
In this section mathematical models are presented for 
the failure rate or hazard function. Recall that the hazard 


may be defined as follows. 


Definition. Suppose that the time to failure, X, is a 


random variable with distribution function F(x), where 


F(0) = 0; the latter possesses the density function f(x), 
mix) = 4Е/ах, such that for any positive x, 
x 
sí e MEUM SA 50 
0 


Then the hazard function, or failure rate at age X, is 
given by 
£ (x) 


h (x) = I- FGE) ° (2.2) 


The interpretation of h(x)dx is that it is the con- 
s TOTI probrbility of failure in the interval (x, x+ dx), 
25: ılo nas been no failure up to age x. 


Express the hazard as 


:—-—- 
MS re 
h(x) dx = fry = - a(los[1 - F(x)]} , 


it then follows after integration that 
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x 
F(x) = 1 - exp[- f h(y)dy] . (2.3) 
0 


Thus if the hazard is specified, so is the distribution 
function, and conversely. 


Note that if 


h(x) 


li 
> 
гё 
о 
о 

[ ^^ 
Ke 


then 
F(x) 


- е ; OR, (254) 


SO a constant hazard function implies the exponential dis- 
tribution of the random variable X, and conversely. 
Obviously a constant hazard representation does not 
describe the bath tub hazard shape of Figure l, nor does it 
represent a situation in which hazards decline, possibly 
because design defects or "bugs" are occasionally removed. 
Here are two hazard representations likely to be useful for 


such purposes. 


1. A Bath Tub Model 


Define the random variable Z in terms of X, X 


being exponentially distributed with mean m. as follovs: 


2 = G(X) = XL(X)R(X) 
or 2) 
= Xp(X), MIK) = L(X)R(X) 
where 
а ОШ ::“ concave in x, L(0) € 1, L(əe) € 1, 


50255 is Convex in X, R(0) = 1, R(0) > R(o). 


r 


oo As 





m 


Then the hazard of Z may be made to exhibit a bath tub 
shape, as in Figure 1, by proper choice of the functions 1 


апа В. 


Example. Suppose 


L(x) - ER 5279000 x 
1+ ах — (2.6) 
= 1 


Clearly, 


ax E ох X 


qa ux = Tg IE BE Dr ox гҹ Ex 


is a monotonically increasing function of x. Furthermore, 
choose a large (e.g. a = 10) and 8 small (e.g. 8 = w 2 
Then it is intuitively clear that (i) small x-values 
transform into even smaller z-values, e.g. x = l corre- 
sponds to Z = 0.91 and x = 2 corresponds to z = 1.9, 
but (ii) this effect dwindles as x increases, so x - l0 
corresponds to z = 9.8 and x = 50 to > = 47.5 and the 
z-values closely resemble the x's percentage-wise, but 
(iii) as x increases still further the z's do not follow 
Saat: x = d corresponds to 2 = 500. This suggests 
that if x is a value assumed by X, that Z shares the 
properties of X in mid-range, i.e. for intermediate 
x-values, but differs from X by having a disproportionate 
probability of assuming small values (near zero), or large 


values (near, but less than, 1/8). Thus the hazard of 


l4 





Zell appear"to be a "bath tubbed” version of X, 
particularly if X is exponential. 

We focus attention on the representation (2.6) in what 
follows, mainly for analytical and computational convenience. 
Se coumse there are many other possibilities, such as 

L(x) = 1 - e 

(27) 

R(x) = PX 

these latter may be adjusted to provide sharper-edged tubs 
than can (2.6), but iteration of (2.6) may be induced to 


accomplish the same purpose. 


2. A Decreasing Failure Rate Model 


Define the random variable W in terms of X, X 


again being exponential with parameter joe: 
W = XT (X) (2.8) 
where T(x) is an increasing function of Xx, L(0) = 1. 


Then the hazard may be made to exhibit a decreasing 


behavior. 


Example. Suppose 


ui EX, CR D NON +. (2.9) 
Then 


2 = EX) (2.10) 


is monotonic, and small x-values lead to comparable z-values 
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y 


(especially when c is small), but larger x-values are 
"amplified" by l + cx to yield increasingly large z-values. 
Attentron will be focused upon (2.9), although other 
possibilities exist that accomplish the same purpose, namely 
that of lengthening the right tail of the distribution of 
X (simulating outliers, for instance) while leaving the 


body of the distribution virtually unchanged. 


B. MATHEMATICAL PROPERTIES OF THE "BATH TUB" HAZARD MODEL 
Various analytical properties of the previously described 

models will now be recorded. These provide useful insights 

into the behavior of the random variables Z and the under- 


lying (generating) variables xX. 


l.  Monotonicity; Quantiles 
It is convenient to focus on Monotonic increasing 


transformations, i.e. if 


z = G(z) = 2012) (210) 


then in order that the above function be monotonically in- 
creasing, dz/dx > 0. Observe that logarithmic differen- 


eration of (2.10) provides 


d27-- $' (x) 


umsthus  dz/dx » 0 if and only if 


o" (x) 
+ $ x) > 0 (2.12) 





x | 
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Alternatively, the condition is, in terms of L(x) and 


R(X); 
T L' (x) R' (x) 


It is easily seen that the important example (2.6), 


n ax : 1. 
2-— r a 1 + Bx ' 


yields a monotonic relationship between z and x. The 
fact that this transformation can be easily and explicitly 
inverted (solved for x win terms of 2) will be exploited 
subsequently. 

Of course if z(x) is monotonically increasing then so 
исми (әли СБ тлуегБе function. The events (Z < z) and 


NA are equivalent, and so 
P(Z € z) — PİX € x(z)) , (2.15) 


from vhich it follovs that if 5. = х(р) 15 Һе р-“ 1002 


auentile of X, 1.0. 


BR =P, (2.15) 
then 


P{Z < z(p)} = P{Z < z(x(p))} = p (216) 


and so 2(р), the p:100% quantile of Z is simply obtained 
from 


zo =p) d(x(p)) € xİp) Lİ(Xx(p)) Rİ(x(p)) (2517) 


In other words we very easily translate from (points on) 
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the inverse distribution of X to the inverse distribution 
of 2. Explicit representation of the distribution of 2 


is however, not often easily possible. 


2. Hazard and Density Function Relationships 


In order to investigate the relationship between 


the hazards of 2 апа X, begin by writing 


“0 
Bern) = 1 - ехр[= f h, (u) du] (2.18) 
0 
or 
x(p) 
h. (u) du = - £n(l-p) 
0 x 


Now differentiate with respect to p to find 


d3055 L 
h. (x (p)) ER - 15 (2.19) 
or 
d ү | 
h (x(p)) = dx (p Tep = f (x(p)) ° I-5 ; (2.20) 


here ns and f. are the hazard and density functions of 
85777 ОК The relationship (2.20) holds for any distri- 
Button, of course. 

Differentiation of (2.5) reveals the connection between 


h and A. Erom 11) 


dz(p) . l $' (x(p)) | dx (p) 
dp cc ‘ats MICE | dp — 


[o(x(p)) + x(p) 9" (x(p)) ] E 
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From (4.11), applied now to the z-hazard, there results 


| _ 1 | 
Rte) © hp) [$04(0)) *x(9) 0" 0:(р))) , (2.22) 
SO 
мз. l 
h,(z(p)) = ho (x(p)) Moe (223) 


Multiplication of both sides by l-p then shows, in view 


of (4.11), that the density functions are similarly related: 
£ (z(p)) = £_(x(p)) 
E x ф(х(р)) + х(р) ф (х(р)) ^ 


Example. 


X is exponential(\). Then 


_ A 
n, (x (p)]) = 4001 + x(p) p x) (2.24) 


Now use the specific (x) of (2.6): 


= ax l 
OR Er é 


or, in terms of logarithms, 


in ф(х) = 4n ox On (l + ах) - ®п(1 + Bx) ; 
SO 
MESE ^. es ша l — aBx (2.25) 
ф(х) x l + ax 1 + 8x х(1+ ах) (1 + Bx)’ ` 
апа 
ф(х) +охф'(х) = ф(х) p ; (2.26) 
finally 
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2 2 
5:77 7 Oa y) (1 + Bx(p)) 
h, (z (p)) – ax(p) [2 + (G + B) x (p) ] (2.27) 


Although this expression is not quite explicit, qualitative 


properties of h, can be deduced rom it. 


WEDA О, х(р) = - š £n(l-p) + 0, and hence 
Ee (2.28) 
2 205007 ` 
Or 
lim x(p) h,(z(p)) = > (2.29) 
р» 0 


since for p - O, 
z(p) + ax“ (p) 
and hence 


x(p) = [z(p)/a]1/2 (2.30) 
there results 


À 
По (2(р)) EE 
2 572575 a 


Or 





lim Vz(p) h,(z(p)) = — er) 
DELI 2/a 


This shows that h, (z) + o as Z > 0, creating the left- 
hand end of the bath tub of Figure 1. 


Re polo =x(p) to, ando z(p) + 1/6 so 


E x” (p) 


Һ. (2(р)) ~ À Тот: 


Or 
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К ш шы (5) = ag? 257 
o co [x (p) ] 2 7 р = a + B ( e ) 
For p>1 
u ир) 
1 В2(р) aB  x(p) 1727 
so 
а + 8 L 
x (p) sl 8z(p) (2234) 
and thus 
В 1 
h_(z(p)) ~a (=) E. 2 
Z 0 ae az)” 


Once again it appears that the hazard rises rapidly, this 
time as x(p) t^» and z(p) 1 871. the other end of the 


bath tub is thus fashioned. 


FED p= 1 — Be? then x(p) = aet Then 
E A 
h,(z(l-e )) = ——— (2.36) 


a/À [2 * (a + B)/A] 


The bath tub effect is presumably achieved by choosing o 
large and 8 small. Let o +0 and B > 0 independently 
in (2.36); it is clear that the limiting value of the 
hazard is A. This indicates that the hazard is (approxi- 


mately) À for middling values of 2. 


ҺГ iron Formula for a Hazard 


The expression (2.6) leads to the relationship 


2 
E ax (р) 
z (p) ~ (I + ax(p)) (1 + 68x(p)) ^ (2.37) 
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süd one latter may be explicitly inverted by solving a 


quadratic equation. The result is 


V = 
х(р) = 5566 + and —— a(l 8z (p)) (2.38) 


Now a direct differentiation of this expression and invo- 


cation of (2.20) produces the expression 


һ (z) 
7: 2180000200 (a48)-2° + 4az (1-82) 
MA 17756 ‘(a + 8) 


E 15-50) z +4az(l-Bz) 


(a + B) z + daz(l - Bz) 


(2.59) 


This form, while explicit, provides no particularly useful 
insights; the bath tub end shapes already noted in (2.31) 
апа (2.35) сап be deduced directly from (2.39). 

Some graphical plots of h, are presented below. They 
illustrate the behavior of the present hazard representation 
in a more understandable fashion than does the formula 


itself. 


4. AN Explicit Formula for the Failure Time Distribution 


Because x and z are monotonically related 


through (2.37) we have 
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Hİ < z} 





_ (a+B)z + V(a+B)"z” + 4az(1-8Z) 
BBE— 7605 у Tm 


200 = BZ 


Il 
FU 
>< 

| ^^ 
e 


| (a+B)z + V(a+B) z* + 4az(l- 82) 


1 — exp % TTS 


(2.40) 


Again the explicit formula seems unproductive of insights. 


C. MATHEMATICAL PROPERTIES OF THE DECREASING FAILURE 
RATE MODEL 
1. The Hazard Behavior 
The expression (2.23) E applied to deduce the 
hazard function of the representation (2.8), advertised to 


produce a decreasing failure rate. There we specified 
1877-7700 75 1 € cx , (2.41) 


52 “rus, from (2.23) 


h (v(5)) - veə TROIS” Ir) (2-42) 
Qualitative properties follov easily. 
ik = 0, x(p) + 0, and 
hg W(P) ~A < (2453) 
Thus the hazard is approximately À for small z. 


017 İp) ? and 


wip) ~ c{x(p)1- 


SO 
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h, (w(p)) ~ —— , (2.45) 


2Vcw(p) 
which clearly decreases, as claimed. It may be inferred 
that the distribution of W appears nearly exponential, 
but possesses an extraordinarily long right tail--possibly 


the result of outliers. 


2- Explicit Formulas for the Hazard and the 


Distribution Function 


Direct solution of the quadratic equation 


wc) ate) a= X + ex” 
presents 
_ /l + 4cw (p) - 1 
x(p) ш r (2.46) 


which, when differentiated, leads to 


1 h (x) 


h,, (W) ———s ”X 
vl + 4cv 


E 5. ( 2.49) 
1 + 4cw 
TU distribution function is 
enz ai el RR 
(2.48) 


1 - exe | - À GA | 


| 


This distribution bears a close family resemblance to the 
Weibull distribution 1 - F(w) = exp{-k/w}, especially for 


large (right tail) values of w. 
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D. AN ALTERNATIVE "BATHTUB" HAZARD REPRESENTATION 

The simple parametric model (2.6) leading to a bathtub- 
shaped hazard is by no means the only possibility. We next 
describe another simple approach. It is that of defining 
a hazard function having an appropriate shape, and then 
deducing the corresponding distribution function, anda 
procedure for sampling from it, rather than proceeding in 
reverse order, as before. 


Let the hazard be of the form 


lize ce giz) mese” k(z) , (2.49) 


vhere g(z) 2 0 is a decreasing function of Z such that 
lim, , , или ила (2) is an increasing function of 2, 


such that (preferably) k(0) 0 and k(9) = o. Such a 


function can yield a bathtubbed hazard. 


Example. 


h(z) = + BZ + A (2. 50) 





A B, &, à all positive. 
Clearly, (2.50) has a generally "bath tub-like" appear- 


ance, since 


Mie ——— pg (DB i) 
(2 + a) 
if 
= - + B < 0, then h'(0) < 0 25) 
vhile for 
A 
- À _ 2.53 
z > 20 = a ( ) 


HZ) 2⁄0. 
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Detailed behavior is adjustable by choice of the 
parameters. Now the distribution function of time to 


failure, Z, is obtained from (2.50): 











2 
C P{Z > z} = expl- f ках | 
0 
T a 
= exp | - / e z t BX + À) dx | 
= exp |- [ a mna + =) +32 + ali 
A 
_ 4 _ (2.54) 
= Fi(z) F5(z) F4(z) , 
where 
s | a NA 
Hü ngam - 
F, (z) = exp (- > 2?) (22555) 
= TAZ 
F, (2) = е : 


All of the above are recognized as being the complements of 
distribution functions. In effect, the distribution of 2 


is that of the minimum of three independent random variables: 


PÍZ » z) = P{X, > 21*Р1Х, > zI-P{X, 2 ZJ; (2156) 
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X, having the distribution F, - 1- F. (п = DU. 

This fact leads directly to an easy procedure for simulation 
of 2 by simply obtaining the smallest from among the 
realization of Xi, Xə, and X4. The advantage of the pre- 
vious method, based on (2.6) for instance, is that only one 
realization--that of an exponential in that specific example 
--leads to the realization of Z. This is not only computa- 
tionally attractive, but seems to facilitate the application 


of such Monte Carlo variance reduction techniques as control 


and antithetic variables, cf., Hammersley and Handscomb [5]. 
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III. OBTAINING SPECIFIED HAZARD BEHAVIOR 


BY SIMPLE SAMPLING 


The development of the last section illustrates one 
manner in which hazard behavior may be conveniently repre- 
sented and simulated. We now show how such behavior may 
alternatively be obtained by simple simulation, i.e. from 
one realization of a basic (possibly exponential) random 
variable. 


Refer to (2.5), in which 


Z = G(X) (3:1) 


and, if G(-) is monotonically increasing, 


z(p) = G(x(p)) , . 665 


4620 x(p) being the p-100% percentiles of 2 апа X, 
respectively. Then the counterpart to (2.23) that' results 


from differentiation of (3.2) is the expression 


T 
һ„(2(р)) = hy(x(p)) == = hy (x(p)) Tdz/ax) (3.3) 


Consequently, if one specifies h, (2) as a suitable func- 
tion of the "time" Z, and specifies the distribution o£ 


the stochastic variable  X--and hence its hazard, h. --there 


results a differential equation for z(x) = G(x): 
dz | 
n, (z) = h. (x) ; (3.4) 


integration then provides the desired transformation, G. In 


other words, we seek 2(х) satisfying 
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ә х 
f h (u)du = / h (v)dv , 5; 
0 De = 


which can sometimes be carried out in a useful closed form. 


Example 3.l. Refer to the example of Section II, wherein 


D. is given by the expression (2.50) and we assume that Х 


is exponential, so n. is constant. Then in order to deter- 





mine G(x) = z(x), solve the equation 
ы x 
/ = + Bu + A| du = f dv 
0 олш 
0 
Or 
д %п(1 + 2) + 22° + 12 = x (3.6) 


Closed-form solution of this expression for Z in terms of 
X is of course impossible. One possible approach is purely 
numerical: find an approximate solution, za (x), e.g. the 


appropriate solution of the quadratic 


z^ + AZ = X (3:77 


təliq 


and then correct the result by a fev Nevton-Raphson itera- 


tions. In other words, put 


NT ES >; 


5 (3.8) 


zi (X) = 


now apply Newton to obtain an improved solution 


A 2n(1+2,/0) 
2900 * 74737 a 5*4 S) 
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which will be feasible if 0 < 22. The process can be 
iterated (the numerator will change after the first iter- 
ation). If one wishes to use this model it may actually be 


desirable to start by solving 


20000 à) z - x = 0 (3.10) 


NU 


for 21, in which case the numerator will not be as shown in 


(3.9); convergence may be more rapid. 


Example 3.2. Change the hazard representation of the previous 


example as follows: let 


ШИШ) = — + Bu + à; a 0) 14) 
(у + a) 
then 
h (udu = А2 + B z” Toz 2) 
0 2 a(z + a) 2 | 3 


Nov it is necessary to solve 


Az 


alz + a) Ы 


Кю 


250 . (0-5 
uc. the cubic 


> = + [a = + MS + [Š + аА - x]z - ax = O , (3.14) 


which can be carried out, at least formally, in closed form. 
Once again an iterative solution that begins by dropping the 
20505 term? solvinq the resulting quadratic for Z (x), and 
then continuing along the Newton-Raphson road may be, 
successful. Further investigations of these ideas should 


be conducted. 
30 





IV. COMPUTER SIMULATION AND ESTIMATION PROCEDURE 


A. SIMULATION AND NUMERICAL RESULTS 

In previous chapters an analytical model was described 
for the failure rate function; useful Ed were also 
derived from the model (2.39) and (2.40). Before the model 
is used in realistic situations, it will be convenient to 


build a computer simulation model for model validation. 


jJ. Algorithm 


First a very basic simulation model was built for 
determining the general shape of the failure rate function 
associated with parameters о, Band À. In the simulation 
model, a was selected to be 1.0 and 2.0, B was selected to 
be 0.05, 0.01, 0.005, 0.001 and À was selected to be 0.1. 
These values were picked arbitrarily, it is only stipulated 
that a is always greater than 8. The general algorithm 
of the simulation model is shown in Fig. 2. 

The hazard function is calculated according to the 
model (2.37) and the system logic function (2.40). The 


results were shown in Fig. 3 and Fig. 4. 


2. Some Comments 
These simulation results show that: 
a. Parameter a is effective when z values relatively 
have small values. That is, it influences the early 


failure period. 
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b. Parameter 8 ¡is effective on the relatively bigger z 
values. That is, it describes wearout failures. 

c. Parameter À has little effect on the shape of the 
curve; it is a scale factor. 

d. The last important observation is that the z values 


are limited by the parameter 8 such that: 
1 
< = 
2-8 
When z equals 1/8, h, (z) goes to infinity. 
B.  EMPIRICAL DENSITY FUNCTION AND ESTIMATION PROCEDURE 
In this section, an estimation procedure for parameters 


a, B, À is defined and the procedure of finding the 


empirical density function is described. 


lu Re E op 
Suppose that Му, N, à and N are defined as 


0 
follows: 
N; = the frequency of data points for each time interval 
between a; and n zer eguüuals l, 2/3; eK 
N = the total number of failures where, 
k 

DT d. E 
A = the length of each interval 

А = b; E. 


sis 





No = the total number of survivors by the beginning of 


each interval which can be shown as follows: 


N ie ael 
N = 
° к= 
N- ) N, 000000 k 
, j=1 J 
where E. - is the number of failures before 


i=l 
interval i, or N - LA N. is the number of 


failures after interval 1. 


IE h_(z) is defined as the density function of 
hazard, then the relationship between the frequences of the 
failure data and density function of hazard can be given 


as follows: 


De 
I 
N.N) fun (z)az m 
a. 
SÉ 
Шеге j = l,2,3,...,k and у = (1-1)А, b; =. JA ee has 


point some approximation can be made in expression (4.1), 


such that 
Ba 
n az Ки ии ann (4.2) 


II 
a. 
l 


for small valves of А = b-a. Then this approximation (4.2) 
is substituted in the expression (4.1), it turns out as 


follows: 


Noe oN 


: ) 
m 0 h, (14) (4.3) 


oF 
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=> 


(1А) = 





Ак (4.4) 


for the general case, it will be: 


7 
N Tore 4 
h (id) = (4.5) 
ae ка 
1-1 “A 
(N- ) N.) S EE UE EON ики 
-. 


The expression (4.5) vill be used in the calculation of 
an empirical hazard function and also used in the proposed 


estimation procedure. 


2. Estimation Procedure for Model Parameters 
Tvo approaches can be used for this problem. The 
first idea is to use the relationship between a sample pth 
percentile and the related probability of the pth percentile. 
A subsequent idea is to approach the problem as a nonlinear 
least square estimation problem for a, 8, A: pick a, 8 
À so that the values obtained minimize the sum of squared 


errors in an objective function. 


ag boot ley Approach 
Suppose PY (x) is the cumulative probability 


Zimekion of the exponential distribution. The pth percentile 


x(p) is equal to the value such that: 


p = Р„(х(р)) = 1-0 MP (4.6) 


BE 





Тһеп 


ЧЕ) = UE (4.7) 


If z(p) is defined to be the pth percentile 
25 bachtub model (2.37) then z(p) and x(p) have a 


definite relationship with expression such that: 


ax” (p) 
- (1 + ax(p)) I + Bx(p)) 


Ey substituting (4.7), 


aí- + 2n(1-p) }* 





ZE) (4.8) 


E we E Lanto) || д + 8 [es in a-p | 


Or 
2“. 
(a/A )€ (p) 
2(р) = ШОК туткыр (4.9) 
{1 + t є(р)} 1 + Š e (p)] 
where e(p) = - £n(l-p). 


In the last equation (4.9), there are three un- 
known variables which are a, B, A. If three independent 
equations associated with expression (4.9) are defined, 
Clearly, they can be solved for the unknown parameters. 
Actually, this can be easily done, for using different values 


Of percentiles such that: 
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(a/3*) * (b) 





2(Р.) = 
L a 8 
(1 + Zeil {1 t с(Р1)) 
| (о/л2)є“ (р„) 
z (P EE 
2 11 A NE ы, 
(4/12) (p,) 


{1 +2 e(P,)| (1 + È e(P3)) 


vhere z (P1) and € (P1) are the known values associated 
with Р.. Actually, to get the percentile values, two kinds 
of approaches can be made. First they can be computed 
directly from data using the simple statistics method, such 
аз 2/3) is (approximately) z(i/(n+1)). Second they can 

be computed from the empirical density function. Generally 
Ehe choice depends on the form of the data that are avail- 
able. 

To decide for the effectiveness of this type of 
estimation, another basic computer simulation is made by 
modifying the first computer simulation algorithm. This 
algorithm is shown in Fig. 5. 

The simulation was run four hundred times and 
in each run the sample size was assumed to be n = 50. 
Initially the parameters a, 8 and ^ were taken to be, 
551000 1.0, 0.05 and 0.1. The percentiles Pi, pə 


and P were used in each run such that 0.1, 0.5 and 0.9 


3 
respectively. The estimation results are shown in Table I. 
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As can be seen, the estimated parameters a, 8, 
À 2: “que different values in all runs. To judge their 
sampling distributions, histograms were drawn separately; 
Figures 6, 7 and 8 are the histograms of the estimated 
values of a, B, А, respectively. 

It is suggested by Figures 6 and 7 for A and 8, 
that normality may be assumed with sample means 0.9246, 
0.05488 and sample variances 0.001065, 0.0002246, respec- 
tively. Their mean, median, trimean and midmean are pretty 
close to each other. But, in the case of a, the histogram 
(Fig. 8) is a quite different picture which looks something 
like an exponential distribution instead of normal distri- 
bution. The reason is clear, because the negative values 
of estimated a were not taken into account. They are 


eivysieally infeasible. 


b. Nonlinear Least Squares Approach 
The least squares criterion used here can be 


stated formally and generally as follows: 


minimize (Y. - Y.) (4.10) 


1012 
> 
tə 


where N is the number of observations, Y: is the fitted 
value of Y.. In this case, expression (4.10) can be re- 


written as follows 


N ` à 
IMS — (a BAN) 
дел = S 


(1611) 
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TABLE I 


COMPUTER SIMULATION RESULTS FOR ESTIMATION PROCEDURE 


Sample size n = 50; Number of Total Runs m= 400; 


Value of parameters a, B, À Introduced: 


Number of 
Runs A o 


a=1, 8=0.05, A=0.1. 


1 0.0858 2.2255 0.0549 
2 0.06591 1.2745 0.0779 
3 0.08381 9.498 0.0688 
4 0.0881 0.9304 0.039 
5 0.1286 3.7429 0.0498 
6 0.0883 214 0.069 
7 0.0785 2.8281 0.0593 
8 0.0945 0.4073 0.0584 
9 0.089 0.3833 0.0452 
10 0.1344 1.3351 0.0426 
iki 0.0677 oan 0.0729 
12 0.1304 2.4871 0.0512 
13 Ol ORA 0.3211 0.0623 
14 0.069 0.3547 0.0657 
15 0.0775 0.4813 0.0772 
16 0.1618 2.7416 0.027 
17 0.0827 0.927 0.0705 
18 0.0816 1.1663 0.0836 
19 0.0515 0.3777 0.0775 
20 0.1089 1.6439 0.0522 
21 0.0767 0.8667 0.0605 
22 0.0499 0.3414 0.0748 
23 0.0569 0.2118 0.0608 
24 0.0648 0.2689 0.0715 
25 0:1415 AOE 0.0518 
26 0.1636 1.3635 0.0186 
27 0.0947 6.101 0.0664 
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TABLE eon t . 
Number of 
Runs A Ы B 
28 0.1063 075235 02.0623 
29 020525 2.6602 0.059 
30 0218 093412 0.0303 
əl 0817 0. 24 0:0655 
32 0.0843 0.4377 0. 0'993 
33 0.1445 4.6941 0.0442 
34 020737 0272956 0 0665 
55 0:20:45 02.7372 0.0646 
36 0.1108 1.1026 0.0427 
E 0. op 2 7 7, 1.0991 0.052 
38 0.0514 0.1362 0.062 
39 02.1.1685 0.6794 0.0469 
40 0.0954 3.7845 0.205987 
41 0 05063 0.4123 0.0738 
42 0.0684 072903 0.0574 
43 0512005 2.059 020557 
44 0; 1885 7.1595 0.0226 
45 0.0708 0.465 0.063 
46 0.0872 2.164 0206057 
47 0.0834 03568 020567 
48 0.0429 0027 0.0694 
49 0.0428 0557295 0.0639 
50 0.0834 0.6807 0.0489 
5 0.1141 1259849 0.0566 
52 0 1959 L. 1959 0.0441 
53 0.1104 570692 0.0433 
54 0.1426 0.8845 0.0244 
55 0.1084 095775 020505 
56 0.1048 150509 0.0562 
57 QE o 033197 0.0643 
58 0.0607 0.4564 0.0607 
59 0.0724 05153 (NS ST 
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Number of 


Runs 


60 
61 
62 
63 
64 
SÉ 
66 
67 
68 
69 
70 
78, 
72 
T3 
74 
nə 
76 
n 
78 
79 
80 
81 
82 
83 
84 
85 
86 
87 
88 
89 
90 


TABLETI Cont. 
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0.9584 
O. TEIG 
0. 1718 
0. 5779 
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0.167 
10129 
022157 
0.2457 
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1, 591099 
002252 
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6.1674 
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Number of 
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am 1 
Li 
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0251 
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12754 
1, ОК 
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6.2455 
1.4348 
0.5245 
ToS 
0525 
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1.5795 
125169 
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9.4492 
0589 
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0D 
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17,608 
közə 
0.3394 
1500029 
2.0504 
3.2468 

























































































































45 


В 


0.058 

0.0482 
0.0611 
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The goal is to choose values a, 8 and A which minimize 
the expression (4.11). In the linear case, this can be done 
using elementary calculus by taking the partial derivatives 
with respect to a, 8 and À , setting each of them equal 
to zero and solving the resulting linear equations (Normal 
equations). But the present problem is different because 
the expression (4.11) is highly nonlinear in the parameters, 
and indefinite in terms of convexity. This is why a non- 
linear programming technique is used. 

There are a few general approaches to the solu- 
tion of the nonlinear estimation problem. One of them is 
ene direct Optimization approach. 


Specifically, equation (4.9) is considered as 


follows: 
2 
0105€ (2)) 
(Pl = TPE BIT + pe tp) — 
BE p 3 p 
where 
_i 
Eege X 
s. 
DET x 
ES 
70 


In equation (4.12), z(p) and e(p) are known values, p. 


m 
are unknown variables. If it is rewritten as the sum of 
squared errors the objective function will be: 
2 2 
N | 0,05€^ (P.) | 
I 2 i 
à : (1 + TE Tur | 5. 
дел | — 270 
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Now nonlinear optimization problems can be stated such that 


“= 25.0 : 
2:7: - ome Р I 5 
gm E Bok 1 PE P; ) 


Subject to 


HD cerco perform the optimization, i.e. to 
solve this nonlinear optimization problem, the GRG package 
(Generalized Reduced Gradient) [6] was used; it is very 
convenient for this case. Ten runs were made with five 
different initial points, both with analytically computed 
derivatives and with numerically computed derivatives. The 
first four initial points were selected arbitrarily. The 


last one was picked such that: 





Ee 
27: 
2 
p, E xo (4.14) 
À 
_ 7(50) 
oT әлан 
À 


where Z0 is the smallest number in the sample, 2 (50) 
is the largest number in the sample, and A is the natural 
logarithm of p = 0.5. divided by the median of the sample. 
ШООЛА points are shown in#Table II. 

The optimization results of a, ß, апа р 


associated with Pyr Po and 03 are tabulated in 


Maple III and Table IV, respectively. 
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As it is seen from Tables III and IV, the program 
was terminated at different optimal points. Because the 
objective function is highly nonlinear and apparently in- 
definite (neither convex nor concave), it may sometimes 
have stopped at a local minimum rather than at a global 
optimum point. But if the results are compared to the 
3-percentile results, the nonlinear least square estimates 


show much greater accuracy and seem more consistent. 


TABLE LE 


INITIAL GUESS POINTS 





V According to the eq. (4.14) 
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TABLE IV 
Actual values À = 0.1, à = 1.0, B = 0.05 


ə. A a B 
Peer | [Initial 


m 0271202 0721208 |”“676.3800 0.0400 
ET 
TITL 
IV 
V 
n 
Hr 
PDT 
IV 
V 
I 
TI 
Ter 
TV 
V 
I 


































05177 355.029 0.0401 







AiO 1720311292591 | 968.150 0.0401 

1 0. T168 1506.926 0.0400 

. 

. 
HET 
IT 
IV 
V 

| 


8221,23 
1501 
204.078 
414.903 
ORA 


392.550 
1520.067 
1465.586 
1177.650 


omo Il 
00332 
0.0330 
00.332 
050 1 


















00115: 874,957 0.0585 
Iu 
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701 | Initial 
Runs Point 
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IV 0.143710.1441 217770) 26 “20710.0312 0.05 
V 41510. 1427 1070508 132610 0.05101]0.0316 
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C. TEST PROCEDURE FOR PARAMETERS 

In order to determine whether the model reasonably fits 
the data a simple test procedure is applied. It is to 
exhibit the simple comparison of the actual points z(P) 
and the estimated values 2(Р), using the estimated 
parameters. 


Basically Z(P) can be obtained as follows: 


“55:67 (4.15) 
A 


Then if it is substituted in expression (2.37), z(P) will 
be 
^^2 


20 -—- ax (P) 


— — (4.16) 
11867 (Р) E Rx (PJ) 


In equation (4.16), all variables are known so estimated 
quantiles can be easily obtained. 

In Fig. 9 actual values and estimated values are plotted. 
Also Fig. 10 displays residuals, which are the difference 


between actual and estimated values. 
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V. NUMERICAL APPLICATIONS TO TWO SETS 


OF REAL DATA 


In this section the two failure data sets studied come 
Bream these sources: Oral irrigator (7) and human life (8). 
They are used for numerical examples in that these dáta are 


fitted using the model (2.6). 


A. ORAL IRRIGATOR 

The data used was obtained from the Commun. Statist.- 
"neor. meth., Colvert and Bordman [7]. The data was col- 
lected such that 100 oral irrigators were placed on the test 
meS terminated at 700 time units. During the life test, 98 
oral irrigators failed; another 2 oral irrigators survived. 
The ordered observed times to failure of the oral irrigators 
5 tabulated in Table V. 

Using the data of Table V, the parameters a, 8, and A 
are estimated by means of the 3-percentile approach. To 
demonstrate the differences between estimated values, 
various different combinations of p values were used. The 
results are shown in Table VI. 

Examination of Table VI indicates that the 3-percentile 
approach may produce estimates having great differences for 
different percentile values, except here in the case of the 
first two combinations. Also the first two estimations seem 
to have accepatable limiting age. Moreoever, the last two 


combinations do not fit well, as judged from the residuals. 
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IABLE V 


TIME TO FAILURE OF ORAL IRRIGATORS 


i75 neo2 15.57 
157.39 ии 25905 
32.44 20957 20:33 
20.99 59.08 86.45 
88.50 103.06 120511 

122.28 122.61 142.27 
142.98 148.29 137.933 
ШОО. 72 16979 233.64 
242.07 256.86 283.95 
288.94 295.48 339.46 
502.61 369.47 414.78 
426.71 459.62 468.64 
469.09 476.42 490.06 
29.3.67 494. 38 418.32 
352.29 334.62 37.0 
E395. 22 309.11 647.91 





TABLE VI 
ESTIMATED VALUES OF a, 8, A BY USING 3-PERCENTILE 





0.00153 0.00634 0.001028 97276 


0.00158 0.00710 0.001018 982.31 
0.00239 0.01835 0.000228 4385.96 
0200273 0.04514 0.000077 12987.01 
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Next the nonlinear optimization approach was applied to 
estimate parameters a, 8 and À associated with Pyr Po 
апа EY both using the derivative and uiu derivative by 
taking into account 98 ordered values. The results аге 
tabulated in Tables VII and VIII for Pyr Por P3 and À 
27— respectively. Some initial points in Table II are 
used. Table VII indicates that the estimates obtained from 
the nonlinear estimation method are more consistent than 
those from the 3-percentile approach. Also, the boundary 
points, 1/8, seem reasonable. Finally, comparison of the 
actual values of z(p) and estimated values of Z(p) indi- 
cate that the fitting is reasonable. The residuals of 


fitted model and the comparison of the actual and estimated 


values are plotted in Fig. 11 and Fig. 12. 


В. HUMAN LIFE (MORTALITY) DATA 

The data used in this example was obtained from the 
1969-71 life table [8] for white females in the United States. 
The table has been prepared from a history of 100,000 persons: 
the number of surviving, and the number dying has been given 


for each age interval. Look at Table IX. 


177 mation by Using 3-Percentile Approach 
The life data is first fitted to the model (2.37) 


by using the 3-percentile approach. Certain pth percentile 
values are selected and used in the estimation process. 
These results are tabulated in Table X. Investigation of 


the last eight combinations of percentiles in these tables 
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Endycates that the limiting age, 1/8, is approximately 110, 
which is reasonable. However, the first combination of 
percentiles gives a good fit in terms of the sum of the 
squared errors. But for these limiting age is reduced to 
about 95 years. This apparently means that the fit of the 
model does not represent the age above 95. Thus it can 

be said that there is a trade-off between 1/8 and the sum 
of squared errors. The present model simply does not seem 


to fit the mortality data very well. 


2. Estimation By Using Nonlinear Least Squares, and 

Application of a Nonlinear Programming Algorithm 

In the previous section we discussed the fact that 
one of the problems encountered in the 3-percentile estima- 
tion procedure was the high uncertainty of fitting, as 
measure by the sum of squared errors. In order to reduce 
this variability, a nonlinear estimation process is applied. 
Previously it was noted that the GRG package can be used 
with analytically computed derivatives and without 
derivatives. If it is used with analytically computed 
derivatives, it is necessary to provide another subprogram 
by the user. Otherwise the GRG package will provide the 
derivative, computed numerically and automatically. 

The derivative of the objective function (4.13) is 


the normal equations such that 


70 





ар 
кшз, 
MES 
qu 
o 
ит 
5 
m 
=> 





TABLE IX 
MORTALITY DATA FOR WHITE FEMALES, 1969-71 


100,000 


99,468 28-29 97,095 73 
98,368 29-30 97 022 
98,303 30-31 96,945 
98,249 31-32 96,864 
98,203 Ss 96,777 
98,164 33-34 96,684 
98,129 34-35 96,583 
98,097 35-36 96,474 
98,068 36-37 96,356 
98,042 37-38 96,228 
98,018 38-39 96,087 
97,994 39-40 95,932 
97,969 40-41 95,762 
97,939 41-42 95,577 
97,902 42-43 95,376 
97,857 43-44 95,156 
97,803 44-45 94,914 
97,743 45-46 94,649 
97,681 46-47 94,358 
97,618 47-48 94,041 
97,555 48-49 93,697 
97,492 49-50 93,325 
97,429 50-51 92,924 
97,365 51-52 92,491 
97,299 52-53 92,022 
97 233 53-54 91,516 





Val 


2 


me + Of # ОЁ Time # Of # Of 
Interval Survivors Dying Interval Survivors Dying 
70 


99,165 








Time 
Interval 


54-55 
55-56 
26 57 
2/56 
58-59 
27-60 
50-61 
pil 62 
42-63 
59-64 
64-65 
55-06 
50207 
67-68 
68-69 
59-70 
70-71 
5-7 
5-/ 
2/4 
14-75 
75-76 
FO 77 
77—78 
78—79 
50 
80-81 
81-82 


# of 
Survivors 





90,970 
90,383 
89,751 
89,071 
88,341 
87,560 
86,726 
85,835 
84,882 
83,860 
82,762 
81,579 
80,304 
78,929 
77,443 
75,836 
74,101 
72,239 
70,246 
68,105 
65,793 
63,290 
60,597 
57,725 
54,686 
51,499 
48,182 
44,748 





















TABLE IX Cont. 


155: Time 


Dying Interval | Survivors 








0—- 
83- 84 
84- 85 
DOS 
86- 87 
60 88 
89 
90 
91 
92 
ə 
94 


















































89= 
OL 
91. 
2 
3 
94- 95 
——: 
20277 
77 Re 
207 
33-00 
100-101 
101-102 
102-103 
103-104 
104-105 
105-106 
106-107 
107-106 
1.098102 
IO 
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# ОЁ 













41,215 
32,629 
34,040 
30,490 
26,995 
23,570 
20,284 
EE 
14,406 
11,875 
9,613 
7,631 
5,937 
4,526 
3,381 
2,476 
1,780 
1,256 
872 
595 
400 
265 












# ОЁ 
Dying 


3,586 
3,589 
3,550 
3,495 
3,425 
3,286 
3,072 
2,806 
2,531 
2762 
1,982 
1,694 
b 
1,145 
905 
696 
524 
384 
Daa 
195 
135 
92 






















re 
712 
71 
45 
28 
17 
























öl 
il 
26 
17 
Il 

6 





20 
i 7 
16 
17 
18 
19 
Ly 
18 
19 
20 


ESTIMATED VALUES OF 


56 
56 
56 
56 
56 
56 
56 
56 
56 
56 


22 

92 

92 
109 
109 
109 
110 
110 
110 
110 


QOO (9 so So Caen 


3-PERCENTILE APPROACH 


Percentile 
Values 


‚000674 
‚000550 
‚000508 
‚000440 
‚000862 
‚000886 
„000908 
500012 70 
5000 
„000920 
2000041 


ÖL, 


TABLE X 


В, 


O O © © O GIO O ©7242 © 


n 


À 


FOR MORTALITY DATA USING 


1.190 
.0404 
0302 
-0199 
. 16 76 
289 
5500 
-1805 
201) 
‚4188 
7700 


CC C © © OO су оо ое 


2010550 
‚010520 
.010550 
-010550 
-009069 
.009068 
2007067 
200621 
. 008989 
. 008988 
-008988 


5: 
777 
94. 
94, 
110. 
1497 
110. 
VLI: 
111. 
WLL; 
lil: 


zl 


N (ә Юю 6) 6) МЮ МЮ 4 





2 
BE. J: | 
001 151 (1 + o,€(P.)) 1 + 03£ (P.)) | 


2 
0105€ (Ру) 


x z (P) = (17 9905901 793657) (4:17) 


2 
as _ _„ Y fo C ese (P0) - e(P,) P 
902 Ea Has 5€ (P.))^ (1 + p4E (P7) 


2 
0105€ (Р.) 


— coU REIS pre 


2 

EN =, ; | 0105€ " (P1) | 
 ——ә— әм. — | re) r 
оз 1=1 | а * o5e (P) ^ (1:4 36 (520) 


` 


010 c^ (P) 
re) ил 

М ls 110 in the case of human mortality data. 

The results of estimation were shown in Tables XI 
and XII for 0; and a, B, À respectively. Also Fig. 13 
and Fig. 14 demonstrate the actual and estimated values and 
£ 0216 Of fitting. 

Examination of Table XI and Table XII indicates 
that the objective function values (the sum of the squared 
errors) of the fitting are smaller than the objection func- 


tion values of 3-percentile approach. Also the estimated 


с 





^ 


parameter values B and A associated with and 


1 ¡dE 
are very close to each other. But estimated values of q 
shov considerable difference for the various initial values. 
However, there are two important facts to notice. First, 
when Q increases significantly, the objective function 
values remain almost constant. That is, it is not effected 
significantly on the objective function values in this case. 
Second, this analytical model (2.6) does not represent 
deaths beyond the age 95. However, modifications can be 


made ın the model (2.6) for this kind of difficulty. Some 


ideas will be discussed in a later section. 


3. Model Modifications 

In the previous section it was noted that the model 
(2.6) has relatively great errors and does not accurately 
describe probability of death at age greater than 95 in the 
human life example. This means that the hazard function 
h, (2) increases too rapidly in the wearout period. If it 
is possible to slow the rate of increase of hazard perhaps a 
better result can be obtained. 

In Section II, it was stated that the function 


R(x) describes the wearout period in bath tub type curves: 


Get 820 


If a new parameter, y, is defined which is between 
0 апа 1 апа R'(x) is now defined to be the m power of 
R(X): 


uo 





525500 - — 1 2:05: (4.20) 


(1 + 8x)! 





then R'(x) will provide for a slowly increasing, rather 
than rapidly increasing, in wearout period. With this 


revision the model (2.6) now becomes 


N 
il 


G(x) = xL(x)*R' (x) (4.21) 


Or 


= 
ns, o e (4.22) 


(1 + ox) (1 + 8x) 
where a > 0, B > О and 0 < Y < 1. 

However, after we put another variable in the model, 
125771) be too difficult to handle in the previous estimation 
technique for obtaining values of the parameters a, B, A 
and ү because of nonlinearity and indefiniteness of the 
expression (4.22). 

For this reason, ү will be assumed constant in the 
previous estimation procedure. Then it will be computa- 
tionally convenient. Actually, the estimation procedure 
does not require any change. The nonlinear least square 
estimation approach will simply be used for different values 
2” aDble XIII and Fig. 15 demonstrate the parameters 
estimated and the resulting objective function values. From 
Table XIII, it can be easily seen that the new result has a 
smaller sum of squared errors. The best value of Y seems 


to Be near 0.95. 
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TABLE ХЕТТ 


ESTIMATION VALUES OF 


USING CONSTANT ү 


c, 


B, 


À 


Run Value 


Јл 


ON 


~] 


7: 
495 
2 
290 


ə 


‚000834 
0.000512 
0.000363 
0.000740 
0.000685 
0.000592 


0.000319 


due 
S 


з 
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.42950 


.04051 
801928 


206374 


81283 
74169 


17467 


01029 


01162 
‚01314 
201352 
-01554 
202051 


„04580 












4,335. 
5119: 
2,494. 
225 
3,00 
5,029: 


157510: 
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Simple analytical hazard models have been developed 
Ende LEted to situations (data) that exhibit bath tub 
shaped hazard functions. That is, failure rates may be 
Rich at early ages ("infant mortality"), constant at 
intermediate ages, and high again for later ages ("wearout"). 
The procedure emphasizes representations of the inverse 
2575 bution function; simulation is thus facilitated. 

The failure time distributions so derived should be 
useful in analyzing maintenance and replacement policies. 

A least squares technique for fitting the hazard models 


to data are suggested and applied. 
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